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Abstract 


NOMENCLATURE 

The transonic flutter dip phenomena on thin 

ah 

distance of elastic axis from the 

airfoils, which are employed for propfan blades, 


origin 

Is investigated using an Integrated Euler/ 
Navier-Stokes code and a two degrees of freedom 

a oo 

speed of sound 

typical section structural model. As a part of 
the code validation, the flutter characteristics 

b 

semi chord 

of the NACA 64A010 airfoil are also Investigated. 
In addition, the effects of artificial dissipation 

Ch 

damping coefficient In plunging 

models, rotational flow, Initial conditions, mean 
angle of attack, viscosity, airfoil thickness, and 

C P 

pressure coefficient 

shape on flutter are investigated. 

Cot 

damping coefficient In pitching 

The present results obtained with a Euler 
code for the NACA 64A010 airfoil are in reasonable 

c 

chord 

agreement with published results obtained by using 
transonic small disturbance and Euler codes. The 

C 1 

lift coefficient 

two artificial dissipation models, one based on the 
local pressure gradient scaled by a common factor 

c m 

moment coefficient about elastic axis 

and the other based on the local pressure gradient 
scaled by a spectral radius, predicted the same 

e 

total energy of the fluid per unit 

flutter speeds except in the recovery region for 


vol ume 

the case studied. 

h 

plunging (bending) displacement, posi- 

The effects of rotational flow, Initial condi- 


tive downward 

tions, mean angle of attack, and viscosity for the 
Reynold's number studied seem to be negligible or 

la 

polar mass moment of Inertia 

small on the minima of the flutter dip. However, 
they have significant effect on the flutter bound- 

J 

Jacobian of transformation 

ary away from the dip. 

<h 

plunging spring constant 

The flutter dip shifts towards higher Mach 
number as the thickness decreases for symmetrical 

K a 

pitching spring constant 

airfoils, other parameters being the same. This is 
In direct relation to the location and strength of 

M 

Mach number 

the shock. The flutter boundary for a thin cam- 
bered airfoil (propfan airfoil) showed that the 

m 

mass per unit length 

effect of camber is nullified by the effect due to 
reduction In thickness and showed a relatively low 

Qh 

generalized force in plunging 

transonic dip. 

Qa 

generalized force In pitching 

The flutter boundary of a simulated SR5 prop- 
fan typical section model showed a very low tran- 

R e 

Reynold's number based on chord 

sonic flutter dip. However, further studies with 
varying mean angles of attack and mass ratio are 

r a 

radius of gyration about elastic axis 

required to better understand the transonic flutter 
dip phenomena of highly swept propfans. 

r cg 

radius of gyration about center of 

*NASA Resident Research Associate. 

S 

gravi ty 

static unbalance, mbx a 
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t time 

t. T nondl mens lonal time, t a /c 

oo 

u, v Cartesian velocities normalized by the 

speed of sound 

V resultant velocity 

Vp flutter velocity 

x distance measured from leading edge to 

trailing edge, normalized by chord 

x,y Cartesian coordinate system 

x a distance between elastic axis and cen- 

ter of gravity 

a pitching (torsion) displacement, posi- 

tive nose up 

a a angle of attack, deg 

n normal direction of transformed coordi- 

nate system 

\i mass ratio 

Ch critical damping coefficient In 

plunging 

C a critical damping coefficient In 

pitching 

£ chordwlse direction of transformed 

coordinate system 

p air density 

u>h uncoupled plunging frequency, A^/m 

u a uncoupled pitching frequency, 7k a /I a 

Superscripts 

( )' d( )/dt 

< )• d< ) df 

n-1 ,n,n+l time levels 

Introduction 

The requirements of high aerodynamic effi- 
ciency and low noise In the operating range of 
transonic flow has resulted In thin (2 to 4 percent 
thick), highly swept and twisted blades for prop- 
fans. These blades operate at moderate to large 
mean angles of attack. In wind tunnel tests of a 
ten-bladed, highly swept SR5 propfan model, 1 con- 
ducted at NASA Lewis Research Center, the blades 
have fluttered at transonic tip Mach numbers at sea 
level conditions. The aeroelastlc characteristics 
of this model were Investigated in Refs. 1 to 3 by 
using a two-dimensional linear subsonic unsteady 
cascade aerodynamic theory with a correction for 
blade sweep. The correlation between theory and 
experiment varied from poor to good, depending on 
the test parameters. To further understand the 
physics of the flutter phenomena and to validate 
the recently developed linear cascade aeroelastlc 
models for three-dimensional subsonic flow condi- 


tions, correlative studies were continued in 
Refs. 3 to 6 which showed good agreement between 
theory and experiment. However, because of the 
limitations of the employed linear subsonic 
unsteady cascade aerodynamic models, the Investi- 
gations in Refs. 1 to 6 were unable to throw any 
light on the effect of nonlinear transonic flow 
on flutter characteristics of swept, thin blades 
in general, and of the SR5 propfan model in 
particular. 

The present investigation was initiated to 
help address the nonlinear effects on the flutter 
of a propfan in the transonic flow regime by using 
a step-by-step approach. In the first step, the 
investigation is restricted to an Isolated airfoil 
with two degrees of freedom. The isolated airfoil 
section Is selected as the section at 75 percent 
span of a thin swept propfan blade. The unsteady 
aerodynamic loads on the Isolated airfoil are cal- 
culated by using a two-dimensional Euler/Navler- 
Stokes solver. This aerodynamic model is selected 
since adequate three-dimensional, unsteady, cas- 
cade, transonic aerodynamic models for rotating 
propfan blades are not available. Even though the 
propfan blades are In a cascade, which has a desta- 
bilizing effect on flutter characteristics, the 
cascade effects are neglected in this investiga- 
tion; since the major goal of this paper is to 
understand the physics behind the flutter dip phe- 
nomena of thin, swept Isolated airfoils. As a 
byproduct of the code validation calculations, the 
effects of blade shape, thickness, mean angle of 
attack, initial conditions, rotational flow, and 
viscosity on the transonic flutter characteristics 
of a thick airfoil are also investigated. 

The approach being used for the flutter analy- 
sis Is the simultaneous Integration of structure- 
fluid dynamics equations. The loads determined 
from solving the fluid dynamic equations are used 
as input to the structural dynamic equations. The 
resultant deflections are then used as input to the 
fluid dynamic equations to determine the loads 
again. If the amplitude of the deflections grow in 
time, then flutter Is said to have occurred. In 
all methods of flutter analyses, the difference 
lies mainly In the prediction of aerodynamic loads. 

Early published research on the effects of 
transonic flow on propfan blade, or propfan airfoil 
flutter characteristics, is nonexistent. However, 
considerable research 7 ” 24 being conducted on 
flutter characteristics of airfoils and aircraft 
wings In transonic flow. Published results 2 ^” 24 
show that there is a considerable amount of reduc- 
tion in flutter speed In the transonic flow regime 
for highly swept blades when compared to the flut- 
ter speed predicted by linear theories. This 
reduction has been called transonic flutter dip. 
This transonic flutter dip has been investigated 
by uslnq transonic small disturbance (TSD) 
theory;™- 20 * 21 modified strip analysis method; 22 
the Euler equations; 23 and the Navier-Stokes equa- 
tions. 24 Many of these methods, along with flutter 
results, are reviewed in Ref. 25. 

The small disturbance theory and full poten- 
tial flow theories fail to take Into account the 
rotational ity of the flow, which could be signifi- 
cant for regions where curved shocks exist. Fur- 
thermore, potential flow theory may give multiple 
steady-state solutions at a given flight condition 
if the free stream Mach number is sufficiently 
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high. Even though the flutter boundaries have 
been predicted using these methods, the question 
has always been asked, how strong are the effects 
of rotational flow, shock location, and shock 
strength on the prediction of flutter boundaries? 
These effects can only be studied by using either 
the Euler equations or the Navler-Stokes equations. 

$ Flutter calculations using Euler equations 

were presented in Ref. 23. The study showed sig- 
nificant differences In the flutter boundaries, 
compared with those obtained using TSD theories, 

f when the steady shocks are strong and at or near 

the trailing edge. Reference 16 presented steady 
and unsteady pressure calculations using Euler 
equations, and concluded that "the Interaction of 
shock with the boundary layer is strong enough to 
alter the airfoil pressure distribution substan- 
tially." The interaction of a shock and boundary 
layer may lead to boundary layer separation result- 
ing In a flow condition which is totally beyond the 
capability of Euler equations. Also, the Euler 
equations do not show the effect of viscosity on 
the flutter boundary. This means that for most 
accurate prediction of the aerodynamic loads and 
the flutter boundary, solution of the Navler-Stokes 
equations Is needed. It should also be noted here 
that the calculations done so far In the published 
literature have been with 6 to 12 percent thick 
airfoils. With the advent of propfans, consisting 
of 2 to 4 percent thick airfoils, an Investigation 
is required for thinner airfoils. 

In the present paper, the Euler/Navler-Stokes 
flow solver, which is described In Ref. 26, and 
coupled to a binary structural dynamic model In 
Ref. 24, will be applied first to study the effects 
of rotational flow, Initial conditions, mean angle 
of attack, viscosity, blade shape, and thickness 
on the transonic flutter characteristics of a 
swept NACA 64A010 airfoil with the structural 
parameters given In Ref. 20. Next, the coupled 
binary structural dynamic model flow solver will 
be applied to a thin, swept SR5 propfan blade 
section to study Its transonic flutter characteris- 
tics. It is to be mentioned here that the artifi- 
cial dissipation model used in Ref. 24 was based on 
the local pressure gradient scaled by a constant 
factor. Recently, efforts are being made 27 to 
incorporate a second dissipation model In this 
code. This model Is based on the local pressure 
gradient scaled by a spectral radius. Preliminary 
results obtained using this model will also be pre- 
sented for comparison. 

Aeroelastlc Model 


A typical section model of a blade, as shown 
In Fig. 1(a), is used for the flutter study. The 
typical section has two degrees of freedom, plung- 
ing (h) and pitching (a), positive as shown. The 
governing equations for this model are 


mh + Sa + C^h + K^h * 

(1) 

Sh + Ia + Ca + Ka-Q 

a a a a 

where m Is the mass, S Is static unbalance, 

Ch,C a are structural damping parameters, and 
K are the plunging and pitching spring constants, 
0? and 0 are the generalized forces, and ( )' is 
h a _ 2 

d( )/dt. Defining h « h/2b, x a * S/mb, r£ « I a / 


mb 2 , u 2 . K h /m, J . KJl*. C h - 2? h <y,. and 

C a « 2t a w a I a , where b Is the semi -chord, I a Is 
the mass moment of Inertia about the elastic axis, 
(*>h and co a are the uncoupled plunging and pitch- 
ing frequencies, (h anc * Ca are the critical 
damping factors in bending and torsion, Eq. (1), 
can be wri tten as 


h + 


a ♦ 


[ 2t n“h] h 


[# * 5mb 


'x 

•• 

Y 

2 

CM 

$- 

3 


'u r 

2 Q 

a 

h + 

a 

•• 

a a a 


a a 

a 

2 

2 

a + 

2 

a ♦ 

2 

a * x 

4irrtr 


The generalized forces, Oh and Qa are 
related to the lift and moment coefficients, 
obtained from the aerodynamic code as follows: 


Q h - - j pV 2 (2b)c ] 
0 » ? pV 2 < 2 b ) 2 c 

a 2 r m 


(3) 


where c^ is the lift coefficient, and c m is the 
moment coefficient about the axis of rotation 
(elastic axis), V is the resultant velocity and p 
Is the air density. Defining, t * ta®/(2b), and 


V/(bw„), where a® Is the speed of sound and 
substituting Eq. (3) into Eq. (2), the final aeroe- 


lastlc equations can be written as 


h + 


'x 

a 

2 


'x 

a 

H 

wsi 

_ | 

h + 

■ 2M<o h ' 

2 

a + 

<V*u ) 
a J 

(V*u ) 

L a j 


Y 

a 

it 

t r m 

a a 

i 

[Mr 1 
a 

2 

a + 

V* 

a + 

V* 


2c ] M 2 


2c M 
m 


ITJi 


(4) 


where ( )' - d( >/dt, M - V/a<», and p is 
m/irp b 2 . For simplicity, the bar on h and t 
will be dropped for the rest of the paper. 

An implicit time marching technique is used to 
solve Eq. (4) In time. Let n be the time level 
where the solution Is known and n+1 is the time 
level where h and a are sought. Then the first 
and second derivatives of h (with similar expres- 
sions for a) are written as 

h ' . a ., -h > ; h . <n zJ tL - ± ±-J. (5) 

(At) 2 

where At Is the time step In the marching proce- 
dure. Substitution of^these finite difference 
derivatives for h', h", a', a" into Eq. (4) leads 
to a system of simultaneous equations for h and 
a at each time step. 

The lift and moment coefficients (c ]f c m ) in 
Eq. (4), are calculated by integrating the pressure 
distribution over the airfoil surface at each time 
step. The pressure is obtained by solving the 
unsteady, two-dimensional, Reynolds- averaged, com- 
pressible Euler or Navler-Stokes equations on a 
body-fitted coordinate system In strong conserva- 
tion form using an alternate direction implicit 
(ADI) procedure. The formulation has been 
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described in Ref. 26, and only a brief outline is 
given here. All the calculations are performed In 
a body-fitted coordinate system (£, n , t) which Is 
mapped to the Cartesian coordinates (x, y, t). 

Fig. Kb), according to the following one-to-one 
relationship: 


5 - 5<x,y,t> 

n - n<x,y,t) (6) 

T * t 

The Jacobian of the transformation J Is 
given by 

J • 5 x n y ' V* ‘ ( x ? y n - x^) <7) 

and the metrics of the transformation are given by 
the relationship: 

* J * n ; *y " ’ Jx n ; n x " - J *e : n y “ Jx ? (8> 


The quantities R 4 and S 4 represent the dissipa- 
tion of energy due to work done by the fluid in the 
x- and y-directlons respectively. The viscous 
stresses t xx , T xy» and T yy are related to the 
velocity gradients through Stokes hypothesis. 

Since the governing equations, Eq. (9), are 
coupled and highly nonlinear, a stable and effi- 
cient solution procedure is required. In the 
present work, the Beam-Warming algorithm , 28 as 
Implemented by Steger, 2 ^ was followed with modifi- 
cations noted by Sankar and Tang . 26 The viscous 
terms are evaluated explicitly to reduce the number 
of operations needed to solve the governing equa- 
tions. A brief description of the solution scheme 
is given here. 


The governing equations are written at a com- 
putational node ( 1 , 3 ) In the following finite dif- 
ference form: 



At 


+ 



+ 6 G 
n 


n +1 


$J n ♦ 5 s" - c c 0 n (13) 
5 n t 


Once a grid has been constructed the metrics 
of the transformation can be evaluated numerically. 
Standard central differences were used to compute 
the quantities such as x n , y n , etc., and these 
quantities in turn were used In Eqs. (7) and ( 8 ) to 
compute 5 X , £y, etc. At the boundaries, three- 
point one-sidea differences were used to compute 
the metrics. 


In the (£, n, t) coordinate system, the two- 
dimensional unsteady Navler-Stokes equations may be 
written as 


where 


q + F- + G 
£ n 


R c + S 
5 n 


(9) 


q - J” (p, pu, pv, e) 


(10) 


and p is the fluid density; u and v are the 
Cartesian components of fluid velocity; e Is the 
total energy^of the^fluld per unit volume. The 
quantities F, G, R, and $ are given by 

, * *v G ♦ 

h * j 


where 


A ~n +1 -n +1 -n 

Aq^ « q - q 


The operators 6r and 6 n are standard cen- 
tral difference operators. For example, the term 
6 {rF m Is the standard two point central difference 
formula given by ( F 1 + 1 1 - Fj-| j)/ 2 . Djj is an 
artificial dissipation term and’Ts discussed In the 
next section. 


The highly nonlinear terms F and G at the 
time level (n+ 1 ) were expanded by a Taylor series 
about a previous time level n as 


P n+1 . F n ♦ (£)" 

V3q7 


(14) 

G n+I - G n + ^q n+1 

where df/dq and 3G/3q are 4 by 4 matrices^which 
are the Jacoblans^of the flux terms f and G 
with respect to q. 


<n v F + n^G + n q> 

r> * Y T 

G j 

« 


In order to allow large values of the explicit 
dissipation coefficient c £ to be used without 
Instability, and to allow the viscous terms to be 
treated explicitly, the following implicit dissipa- 
tion terms were added to the left side of the dif- 
ference Eq. (13). 


S 


<n x R ♦ n y S) 
J 


-D.r’a., + 5 >J bq n+1 
I SS nn 


(15) 


The terms F and R are the flux and vis- 
cous stress terms along the x-dlrectlon and G 
and S are likewise the flux and viscous terms 
along the y-dlrectlon and are given as follows: 


The coefficient ej was taken to be two to 
three times the explicit dissipation coefficient 
e^. A range of values between 2 and 5 were 
used In the calculations. 


F - (pu, pu 2 + p, puv, u(e + p>) 
G - (pv, puv, pv 2 + p, v(e + p>) 
R - (0, T XX , T X y, R 4 ) 

S - (0, T X y, Tyy, S 4 ) 


The dissipative term, Djj, Is added to remove 
high-frequency errors In the solution at every time 
step, and to avoid odd-even decoupling of the 
( 12 ) numerical solutions. Two models of introducing 

artificial dissipation are available in the present 
code. In model I , 26 the dissipation term is writ- 
ten, in conservative form, as a combination of sec- 
ond and fourth order dissipative terms. This model 


1 


4 



I 


t 


Is based on the local pressure gradient scaled by 
a constant factor. A sensor based on the second 
derivative of pressure, turns the second order dis- 
sipation term In the vicinity of shocks and sup- 
presses the fourth order dissipation terms. This 
sensing is followed only In the streamwlse direc- 
tion. A more detailed description of this model, 
and bench mark calculations are presented In 
Ref. 26. This model, as pointed out In Ref. 30, 
adds positive dissipation in certain regions and 
negative dissipation In other regions depending 
on the flow field gradients In that region, which 
may degrade the solution accuracy in some flow 
problems . 

In model II, 27 the local pressure gradient Is 
scaled with the spectral radius of the flux vector. 
The fourth order difference term Is modified In a 
way that It produces positive dissipative terms. 
This model Is similar to the one presented In 
Refs. 31 and 32 and Is being currently Implemented 
in the present code by Wu. 27 Some preliminary test 
results will be given In the present paper. 

Equation (13) may be written after the addi- 
tion of the artificial Implicit dissipation terms 
given by Eq. (15), In the following operator form 

I + At 6 r — + At 6 — - e r J _1 (6,, + 6 )J 

l 5 aq n 3q 1 « no J 


slip boundary conditions were enforced by setting 
the velocity of the fluid to that of the solid sur- 
face at every time step. At the solid surface, the 
normal derivative of the temperature and pressure 
were also set to zero. The use of the C-grid sys- 
tem Introduces a cut between the airfoil trailing 
edge and the downstream boundary. Along this cut, 
the flow properties were averaged from above and 
below. The body-fitted grids employed in this work 
were such that the far field boundaries were at 
least six chord lengths away. The flow properties 
at the far field boundaries were assumed to be 
undisturbed . 

The Navler-Stokes solver described here may 
also be used to model inviscld rotational flows. 
This requires that the viscous terms In the Navier- 
Stokes equations be suppressed, and that the zero 
tangential velocity boundary condition at the solid 
surface be replaced by a nonzero value, obtained by 
extrapolating the tangential component of the fluid 
velocity from the interior. 

The code Is expected to predict the flow field 
for arbitrary airfoils for moderate angles of 
attack, a Mach number range from 0.1 to 1.2, and a 
Reynolds number range from 100 to 10 7 . This code 
has been used to predict the dynamic stall charac- 
teristics for the NACA 0012 airfoil ,26 transonic 
flow solutions, 36 and stall flutter characteristics 
for the NACA 0012 airfoil . 24 


Aq n+1 - R n (16) 


Results and Discussion 


where 

R n - -At(6 c F + 6 G) n + At(6 r R + 6 S) n - At e F D n 
t n £ n t 

(17) 

The left hand side operator of Eq. (16) was 
approximately factored Into two small operators, 
leading to the following final form: 


This section presents the numerical results 
obtained with the aeroelastlc model described ear- 
lier. First, the code is validated with published 
results. Next, the effects of rotational flow, 
Initial conditions, mean angle of attack, and 
shape and thickness on transonic flutter are Inves- 
tigated. Then transonic flutter calculations are 
done for the SR5 propfan blade typical section 
model . 

Code Validation 


At S . — - e. AtJ' 1 

*• 3q ‘ 


At 6 — - e. AtO' 1 5 

n a5 I nn 


A ^n + 1 n n 

Aq + R 


(18) 


Equation 18 may be solved through the Inver- 
sion of two block trldlagonal matrix equations, one 
corresponding to the ^-direction, and the other 
corresponding to the n-dlrectlon. In order to 
keep the flow solver simple, the boundary condi- 
tions on all the boundaries were explicitly updated 
after the interior points had been updated using 
Eq. (18). 


The Baldwi n-Lomax 33 eddy viscosity model Is 
used to model the turbulent momentum and energy 
transfer. This is a two layer algebraic turbulence 
model with two algebraic equations In each layer. 


Algebraic grid generation routines were 
included by Sankar and Tang26 to make the code more 
portable. Both 0- and C-grlds can be generated. 
However, any grid can be overwritten for the alge- 
braic grid, for example, by a numerical grid. 34 * 35 

The boundary conditions were treated expli- 
citly as follows: at the solid surface, the no 


References 24, 26, and 36 presented selected 
validation cases with the present Euler/Navier- 
Stokes flow solver. These cases represented dif- 
ferent airfoils at different flow conditions. In 
those references, the dissipation model I, which 
was described earlier was used. As a part of the 
validation study of the dissipation model II, 
described earlier, the pressure distribution for 
RAE2822 airfoil Is recalculated. The results are 
compared with those obtained In Ref. 36 with dissi- 
pation model I and the experimental data presented 
therein. Figure 2 shows the steady pressure dis- 
tribution (pressure coefficient, -C p , versus non- 
dimensional distance along the chord, 
x - 0.5*(x/b + 1), see the coordinate system in 
Fig. 1(a)) for the RAE 2822 airfoil at Mach 
number, M « 0.73, angle of attack, a a » 2.79°; and 
Reynold's number, Re * 6.5 million. The numerical 
results with both dissipation models and the exper- 
imental data show very good comparison. These com- 
parisons help to give confidence In the code and In 
the attempt at flutter calculations with the dissi- 
pation model II. 

In order to establish the capability of the 
code for predicting the transonic flutter dip, the 
flutter boundary of the NACA 64A010 airfoil Is 
obtained. The calculations are carried out using 
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the Euler version of the code with both dissipation 
models. A 157 by 40 algebraic C-grld shown In 
Fig. 3 Is used In the calculations. The structural 
parameters are those used In Ref. 20 and are given 
as case A In Table 1. In Ref. 20 t the wing geomet- 
ric sweep is modeled by positioning the elastic 
axis ahead of the leading edge (ah * -2.0, see 
Fig. 1(a)). 

To obtain the flutter boundary, a steady-state 
flow field solution Is first obtained for 0° mean 
angle of attack. The solution of the aeroelastic 
equations (Eq. (4)) Is then started at an assumed 
nondlmens lonal flutter speed, V*(*V/( bco a ) , by pre- 
scribing the initial conditions, h(0>, a(0), h (0), 
and a' (0) . If the amplitude of oscillation of the 
resulting transient solution increases (decreases), 
the value of V* Is reduced (Increased) and the 
transient solution Is again calculated. The criti- 
cal flutter speed (V*p) Is obtained by linearly 
interpolating the logarithmic decrement of the 
amplitude of oscillations for two different but 
close values of V* (one giving oscillations with 
increasing amplitude and the other giving oscilla- 
tions with decreasing amplitude). The flutter 
boundary is presented as a variation of the flut- 
ter speed index (V*p/ >/jr) versus Mach number (M). 

The results obtained for the NACA 64A010 air- 
foil are compared with published resul ts W, ^ *23 
in Fig. 4. Refs. 10, 11, and 21 used transonic 
small disturbance (TSD) theory and Ref. 23 used 
Euler equations in the cal culations . In this com- 
parison, differences due to grids, numerical 
schemes, artificial viscosity, etc., employed in 
the published results are not considered but only 
the flutter boundaries are compared. In the 
present calculations, the Initial conditions, 
a( 0) - 0.1°, with zero h(0), h'(0), and a‘(0) 
are used. The comparison shows that all of the 
codes '0* 1 1 » 21 *23 qualitatively predict the tran- 
sonic dip. 

For Mach numbers between 0.7 and 0.82, the 
flutter speeds predicted by the present Euler code 
are less than the corresponding values from 
Refs. 10, 11, 21, and 23. For Mach numbers between 
0.82 and 0.87 the flutter speeds predicted by all 
the codes are In good agreement. These differences 
can be explained partly by Inspection of the pres- 
sure distribution over the airfoil. For a Mach 
number of 0.8, the pressure distribution (C p ), 
showed a shock appearing near the mid chord. This 
shock travels towards the trailing edge and 
Increases in strength as the Mach number Increases. 
The T SD theories are not able to model the effect 
of rotational flow downstream of the shock which 
could be significant in some cases. Since the 
shock is almost near the mid chord, the flow is 
rotational over a reasonably large surface. This 
effect on a considerably large surface Is totally 
neglected by the TSD codes which could result in a 
higher flutter boundary. Thus the difference 
between the flutter speeds calculated by using the 
Euler code and TSD codes in this region may be 
attributed to rotational flow and entropy effects 
behind shock. However for Mach numbers between 
0.85 and 0.88, the shock Is almost near the trail- 
ing edge reducing the surface over which the flow 
is rotational. Then the calculations with Euler 
code shows same accuracy as of TSD codes. This 
could explain the very good comparison by all the 
codes between M - 0.85 and 0.88. 


However, the transonic recovery portion of the 
flutter boundary predicted by the Euler codes diff- 
ers from those predicted by codes based on TSD 
theory. The recovery predicted by Euler codes 
occurs at a higher Mach number compared to the one 
predicted by TSD codes. Also, the Euler code of 
Ref. 23 predicts the recovery to be at a higher 
Mach number compared to the present Euler code. 
Beyond M - 0.88, the shock strength Is suffici- 
ently high enough to Induce separation. The Euler 
and TSD codes fall to model separated flow behind 
shocks and predict qualitatively the same flutter 
boundary beyond M - 0.88. 

Also, there are significant qualitative diff- 
erences in the recovery region of the flutter 
boundary obtained by the two dissipation models 
available in the present Euler code. The dissipa- 
tion model I predicts a flutter boundary which 
folds over to provide upper boundaries. For a Mach 
number of 0.9, three boundaries are observed. This 
compares qualitatively with the boundaries obtained 
In Refs 10, 11, and 23. Whereas no such folding 
over of the boundary is observed when dissipation 
model II is used. Between Mach numbers of 0.88 and 
0.9 the flutter boundary jumps considerably, and is 
almost the same as the third boundary obtained with 
dissipation model I. This has also been observed 
by Isogal.21 Again, this difference can be 
explained by observing the pressure distribution 
obtained with the two dissipation models. The 
pressure distributions obtained at M = 0.9 with 
the two dissipation models are shown In Fig. 5. 

Both models show considerable differences with 
regard to the shock location and shock strength. 

The dissipation model II predicts the shock to be 
closer to the trailing edge than that predicted by 
dissipation model I. This may be the reason for 
the difference observed in the recovery region with 
these models. 

Figure 6 shows a typical transient response 
solution. Figure 6(a) shows the transient response 
for a point below the flutter speed and Fig. 6(b) 
shows the response for a point above the flutter 
speed. In both the responses, the amplitude of h 
Is larger than that of a. indicating possible 
bending mode instability. 

The above comparisons for the NACA 64A010 help 
to give confidence In the code for investigation of 
the flutter dip phenomena for other airfoils. It 
also shows the Importance of the accurate predic- 
tion of the shock strength and position. 

The code is now applied to study the effects 
of the following: (1) initial conditions, 

(2) mean angle of attack, (3) viscosity, and 
(4) thickness and shape. Dissipation model I is 
used In all the calculations. 

Effect of initial conditions . To study the 
effect of initial conditions, the Initial condition 
on a(0) is increased from 0.1° to 4°. The flutter 
boundary for 4°, together with that for 0.1°, Is 
shown in Fig. 7. For a Mach number of 0.8 the 
flutter speed for 4° decreases by approximately 
30 percent. This has also been observed in 
Ref. 23. However, for Mach numbers of 0.85 and 
0.88 there does not seem to be any significant 
change In flutter speeds (less than 10 percent). 

This again seems to be due to the position of the 
shock. Since the shock is near the trailing edge 
of the airfoil, the effect of rotation does not 
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seem to be of much significance. The results show 
that the Initial conditions coupled with transonic 
flow nonl Inearl ties have significant effects on the 
flutter characteristics In the transonic Mach 
number range, but negligible effect on the minima 
of the transonic dip. 

Effect of mean angle of attack . The effect of 
mean angle of attack was studied by starting the 
• solution of the aeroelastlc equations (Eq. (4)) at 

different steady flow field conditions. In the 
present study, the steady flow field Is calculated 
at two degrees of angle of attack, Instead of 0°, 

» for Mach numbers 0.75, 0.8, 0.85, and 0.88. The 

Initial condition, a(0), Is again 0.1°. These 
results are also Included In Fig. 7. The transonic 
dip moves to left since the shock occurs at lower 
Mach numbers. For Mach number equal to 0.8, the 
reduction In flutter speed Index Is almost 50 per- 
cent. This was also observed by Edwards, et al. 10 
However, the shifting of the flutter boundary does 
not seem to be as much as observed In Ref. 10. 
Again, the effect of the mean angle Is negligible 
at the minima of the dip. 

Effect of viscosity . The flutter boundary cal- 
culations to study the effect of viscosity are very 
expensive. The required computational time Is 
almost an order of magnitude higher than that for 
the Euler case. Hence the effect was studied only 
for a Reynold's number of 12 million. These 
results are also shown In Fig. 7. For this case 
again the flutter boundary near the dip shifts 
down by less than 5 percent. This shift could be 
attributed to change In grid spacing near the air- 
foil since a very fine grid Is required near the 
airfoil for the viscous case. The only significant 
difference was that even with the dissipation 
model I, there was no upper boundary observed. 

The effect of viscosity seems negligible for the 
Reynold's number studied here. However, studies at 
a lower Reynold's number are In progress. 

Effect of shape and thickness . The effect of 
shape and thickness on the flutter boundary Is 
studied next. Four NACA 16-series airfoils with 
the same structural dynamic parameters as those 
used for the NACA 64A010 airfoil are considered for 
the study (see Table 1, case A). This series of 
airfoils Is selected because these airfoils are 
used on current propfan designs. For example, the 
SR5 propfan has 1 6- s e r 1 e s airfoil sections on the 
outer 45 percent of the span. The airfoil sections 
vary In thickness and camber. The section selected 
for the present analysis Is 2.6 percent thick and 
has a design lift coefficient of 0.13. Therefore, 
for the study of shape and thickness, the airfoils 
selected are 16-010, 16-004, 16-(1 . 3 > ( 04) , and 
16-C1 .3X2.6). Here the first two digits specify 
the series, the next digit Indicates the amount of 
camber expressed In terms of the design lift coef- 
ficient in tenths, and the last two digits specify 
the thickness to chord ratio In percent chord. The 
first three airfoils are Included for a step-by- 
step verification of the code and to applications 
for thin airfoils. Again, an algebraic C-grld of 
1 57 by 40 i s used . 

The flutter boundaries for the selected air- 
foils are compared In Fig. 8. Comparison of the 
results for the NACA 16-010 and the NACA 16-004 
airfoils show that the transonic dip shifts to the 
right for the latter airfoil. This Is due to a 
decrease In thickness which delays the formation 


of the shock to higher Mach numbers. Because of 
smaller thickness, the shock strength is also small 
relative to that of NACA 16-010. This results In 
a smaller jump In the flutter boundary In the 
recovery region. For the same Mach number, the 
NACA 16-004 airfoil has a shock appearing behind 
that of the NACA 16-010 airfoil which results in a 
smaller region of rotational flow for the NACA 
16-004 airfoil. The combined effect of the shock 
location and shock strength, and rotational flow 
Is to lower the boundary for the NACA 16-010 as 
compared to that of the NACA 16-004 for the same 
Mach number. The flutter boundaries calculated up 
to M * 0.9 for the NACA 64A010 and the NACA 
16-010 airfoils did not show any difference, indi- 
cating thickness has more effect on the flutter 
boundary than the shape. 

The flutter boundary obtained for the NACA 
16-(1.3)(04) Is also shown In Fig. 8. The pressure 
distribution curves, not shown here, for this air- 
foil showed a shock at a Mach number equal to 0.85, 

which resulted In a drop In the flutter speed index 
compared to that for the NACA 16-004 airfoil. The 

shock Is at the trailing edge beyond M » 0.88, and 

the flutter speed Index starts Increasing again. 

The similarity of the flutter boundary curve with 
that of the NACA 16-004 airfoil Is consistent with 
the observation that the camber effect Is equiva- 
lent to the effect of Initial angle of attack for 
the NACA 1 6-< 1 .3X04) airfoil. 

In Fig. 8. the flutter boundary obtained for 
the NACA 16-0 .3X2.6) airfoil (the propfan air- 
foil) Is also shown. A dip in the flutter boundary 
can be seen near M - 0.925. There after the flut- 
ter speed Index starts Increasing. It is interest- 
ing to note that the flutter boundary between Mach 
numbers 0.85 and 0.95 Is almost similar to that 
obtained for the NACA 16-004 airfoil. This is due 
to the fact that the reduction In the flutter speed 
Index due to camber, like that for the NACA 
16-0. 3X04) airfoil, Is compensated by a decrease 
In thickness. Between Mach numbers 0.75 and 0.85, 
the flutter boundaries for the NACA 16-004 and the 
NACA 16-0.3X2.6) airfoils show the effect due to 
thickness, as was observed for the NACA 16-010 and 
the NACA 16-004 airfoils. A typical response for 
the NACA 16-004 Is shown In Fig. 9. Figure 9(a) 
shows the transient response for a point below the 
flutter speed, and Fig. 9(b) shows the response for 
a point above the flutter speed. In both responses, 
the amplitude of h Is larger than that of a, 
Indicating possible bending mode Instability. 

The comparison of the results from the four 
airfoils shows that (1) for symmetric airfoils, the 
transonic dip moves towards higher Mach numbers as 
thickness to chord ratio decreases, (2) for the 
thin cambered airfoil studied here, the camber 
effects are nullified by the effects due to reduc- 
tion In thickness, and the flutter boundary is 
similar to a thick symmetric airfoil of the same 
series, and (3) the flutter characteristics 
strongly depend on the shock location and shock 
strength. 

After establishing confidence In the code, it 
was applied to Investigate the flutter characteris- 
tics of a simulated SR5 propfan blade. 
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Flutter Characteristics of a 
Simulated SR5 Propfan Blade 


Conclusions 


A typical section model Is constructed which 
simulates the structural properties, mode shapes, 
and coupled frequencies of the SR5 propfan. Linear 
subsonic Isolated airfoil theory^' Is used In simu- 
lating the typical section model. The rotation 
axis (elastic axis) and the uncoupled bending and 
torsional frequencies are chosen such that the mode 
shapes and coupled frequencies are close to those 
calculated for a SR5 propfan blade using NASTRAN 
analysis. The selected properties are shown In 
Table 1, case B. It should be noted here that 
the frequency ratio Is almost equal to one as In 
case A, studied earlier, and the mass ratio Is 
approximately twice that of case A. For the struc- 
tural properties selected, the calculated flutter 
Mach number and flutter reduced frequency are 0.875 
and 0.17. The corresponding experimental 1 values 
are 0.88 and 0.17 respectively. 

The flutter boundary obtained using the 
present Euler code with dissipation model I Is 
shown in Fig. 10. The flutter boundary obtained 
with linear Isolated blade theory Is also shown in 
Fig. 10 along with the velocity Index curve 
(V/(b« a p) versus M) for this typical blade. 

The present Euler code gives a flutter Mach number 
of 0.845. This value Is about 4.5 percent less 
than that obtained from SR5 propfan experiment. 

The prediction of lower flutter Mach number by the 
present Euler code may be due to the following rea- 
sons: (1) the limitation of the typical structural 

dynamic model employed to simulate the SR5 propfan 
and (2) no structural damping Is included in the 
present analysis. 

Figure 10 shows that the flutter speed Index 
drops as the Mach number Increases up to M - 0.9, 
thereafter, it starts recovering, showing a very 
low transonic flutter dip for SR5 propfan blade. 

The Euler code predicts almost the same flutter 
boundary as the linear Isolated blade theory. How- 
ever, a much more pronounced transonic dip Is 
expected when the following aspects are considered. 
The above calculations are made at zero mean angle 
of attack. In a real situation, the propfan Is at 
a mean angle of attack which may cause transonic 
flutter dip. The mass ratio used is 115, which Is 
higher than for the present composite propfans hav- 
ing the mass ratio of the order of 33 to 60. This 
low mass ratio may cause transonic flutter dip 
which needs further Investigation. The cascade 
effects which are found to have destabilizing 
effect on flutter In subsonic flow, may also effect 
the transonic dip. 

The following points support the above points 
when the flutter boundary obtained for the NACA 
1 6- ( 1 .3X2.6) airfoil In Figs. 8 and 10 are com- 
pared. In Fig. 8, a relatively large dip Is seen 
for the same airfoil section. Here the mass ratio 
Is 60, which Is about 48 percent less than that 
used for SR5 propfan blade typical section In 
Fig. 10. Also, the elastic axis position Is one 
semichord length In front of the leading edge com- 
pared to that used for SR5 propfan blade typical 
section. Since here only an attempt Is made to 
simulate a SR5 propfan, these structural properties 
need to be improved for correct prediction of flut- 
ter boundary and correlation with experiment In 
transonic flow for thin highly swept propfans. 


An Euler/Navler-Stokes flow solver Is used to 
obtain the flutter boundaries for a typical section 
structural model with several airfoil sections. 

The application of the code for predicting tran- 
sonic flutter characteristics of thick and thin 
airfoils Is demonstrated. Based on this study, the 
following conclusions are drawn. 

1. The rotational flow effects behind the 
shock have a strong effect on the transonic flutter 
speed depending on the chordwlse location of the 
shock. The neglect of rotational flow effects 
results In predicting a higher flutter speed. 

2. The two dissipation models employed, one 
based on the local pressure gradient scaled by a 
constant factor and the other based on the local 
pressure gradient scaled by spectral radius, pre- 
dicted the same flutter boundary up to the tran- 
sonic dip for a thick airfoil. However, in the 
recovery region, the boundaries predicted by these 
models are different. 

3. The effects of Initial conditions, mean 
angle of attack, and viscosity for the Reynold's 
number studied on the minima of the transonic dip 
seems negligible. However, they have a signifi- 
cant effect away from the dip. 

4. The blade thickness and shape dictate the 
location and strength of the shock, there by 
affecting the flutter boundary. The transonic dip 
shifts to higher Mach numbers for symmetric air- 
foils with decreasing airfoil thickness to chord 
ratio. The flutter boundary, for relatively thin 
cambered airfoils studied here, showed a relatively 
low transonic dip compared to a symmetric thick 
alrfol 1 . 

5. The predicted flutter Mach number for a 
simulated SR5 propfan blade Is about 4.5 percent 
less than that obtained by experiment. This dif- 
ference Is attributed to the simplified aeroelastlc 
model used In the present analysis. The flutter 
boundary showed a very low transonic dip. Further 
studies are needed to investigate the effects of 
mean angle of attack, mass ratio, and elastic axis 
position on the transonic flutter of highly swept 
composite propfans. 
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Table 1. - STRUCTURAL PARAMETERS 
USED IN THE STUDY 


[Distances are nondimensional ized 
with respect to semi-chord (b), 
see Fig. 1(a) for definitions.] 


Parameter 

Case A 

Case B 

a h 

wh/w a 

*a 

r cg 

p 

Ch.Ca 

-2 

1 

1 .8 

0.4898 

60 

0 

-1 

0.928 

0.9644 

0.4769 

115 

0 
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AMPLITUDE 


.00400 
.00250 
.00100 
- . 00050 
-.00200 

.00500 

.00325 

.00150 
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-.00200 




(B) V = 5.0. 


FIGURE 6. - PLUNGING AND PITCHING AMPLITUDES VERSUS 
TIME, NACA 64A010 AIRFOIL, M = 0.88, a h = -2.0, 
V w a = 1-0- x a = 1.8, r a = 1.865, EULER, ALGE- 
BRAIC C-GRID (157x40). 
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FLUTTER SPEED INDEX, V, 



MACH NUMBER, M 


FIGURE 10. - FLUTTER BOUNDARY FOR SR5 PROPFAN SIMULATED 
TYPICAL SECTION MODEL. (CASE B, TABLE 1, DISSIPATION 
MODEL I.) 
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